###################################################################
### Natural Hazards, Social Policy, and Electoral Performance: ####
### Evidence from the 2017 Earthquake in Mexico City           ####             
### Figures and tables in the main text                        ####                  
### This version: February 26, 2023                            ####                   
###################################################################
 
#---------------------------------------------------#
# Preliminaries                                     #
#---------------------------------------------------#

rm(list = ls())

#List of packages for session
.packages = c("tidyverse", "stringi", "stringr",
              "estimatr", "patchwork", "modelsummary",
              "spdep", "sf", "MatchIt", "texreg", "sensemakr")

#Install CRAN packages (if not already installed)
.inst = .packages %in% installed.packages()
if(length(.packages[!.inst]) > 0) install.packages(.packages[!.inst])

#Loading packages into session 
lapply(.packages, require, character.only=TRUE)

# Set Working Directory to wherever files are downloaded
#setwd("")

#---------------------------------------------------#
# Loading Data                                      #
#---------------------------------------------------#

#Main datasets for regression analysis and figures

##Master dataset
master <- read.csv("master_cdmx.csv") %>%
  dplyr::select(-1)

###############################################
############ Main Text -- Figure 6 ############                     
###############################################

#Figure 6.  Relationship between damage associated with the earthquake and the electoral performance of the PRD candidate for governor.

#City Incumbent (PRD)

#Models with distance

model1_gov <- (list(
  "Any Damaged Housing" = lm_robust(prd_change_gov_12_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(prd_change_gov_12_18~distance_damaged_houses_highrisk+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(prd_change_gov_12_18~distance_damaged_multifamily+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('distance_damaged_houses' = 'Distance to Damage',
               'distance_damaged_houses_highrisk' = 'Distance to Damage',
               'distance_damaged_multifamily' = 'Distance to Damage')

plot1_gov <- modelplot(model1_gov, 
                       coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                       coef_map = coefnames) +  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Incumbent (PRD)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + theme_light()

#Models per capita

models2_gov <- (list(
  "Any Damaged Housing" = lm_robust(prd_change_gov_12_18~damaged_houses_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(prd_change_gov_12_18~damaged_houses_highrisk_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(prd_change_gov_12_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('damaged_houses_1k' = 'Damage Per Capita',
               'damaged_houses_highrisk_1k' = 'Damage Per Capita',
               'damaged_multifamily_1k' = 'Damage Per Capita')

plot2_gov <- modelplot(models2_gov, 
                       coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                       coef_map = coefnames) +  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Incumbent (PRD)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + theme_light()


plot1_gov + plot2_gov + plot_layout(nrow = 2)

###############################################
############ Main Text -- Figure 7 ############                     
###############################################

#Figure 7.  Relationship between damage associated with the earthquake and the electoral performance of the PRD and MORENA candidates for mayor. 

#City Incumbent

#Models with distance

models1_mayor <- (list(
  "Any Damaged Housing" = lm_robust(prd_change_mayor_15_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(prd_change_mayor_15_18~distance_damaged_houses_highrisk+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(prd_change_mayor_15_18~distance_damaged_multifamily+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('distance_damaged_houses' = 'Distance to Damage',
               'distance_damaged_houses_highrisk' = 'Distance to Damage',
               'distance_damaged_multifamily' = 'Distance to Damage')

plot1_mayor <- modelplot(models1_mayor, 
                         coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                         coef_map = coefnames) +  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Incumbent (PRD)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + theme_light() + xlim(-2,2)

#Models per capita

models2_mayor <- (list(
  "Any Damaged Housing" = lm_robust(prd_change_mayor_15_18~damaged_houses_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(prd_change_mayor_15_18~damaged_houses_highrisk_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(prd_change_mayor_15_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('damaged_houses_1k' = 'Damage Per Capita',
               'damaged_houses_highrisk_1k' = 'Damage Per Capita',
               'damaged_multifamily_1k' = 'Damage Per Capita')

plot2_mayor <- modelplot(models2_mayor, 
                         coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                         coef_map = coefnames) +  
  labs(title = "City Incumbent (PRD)") +  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Incumbent (PRD)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + 
  theme_light() + xlim(-0.012, 0.012)

#City Challenger

#Models with distance

models3_mayor <- (list(
  "Any Damaged Housing" = lm_robust(morena_change_mayor_15_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(morena_change_mayor_15_18~distance_damaged_houses_highrisk+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(morena_change_mayor_15_18~distance_damaged_multifamily+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('distance_damaged_houses' = 'Distance to Damage',
               'distance_damaged_houses_highrisk' = 'Distance to Damage',
               'distance_damaged_multifamily' = 'Distance to Damage')

plot3_mayor <- modelplot(models3_mayor, 
                         coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                         coef_map = coefnames) +  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Challenger (MORENA)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + theme_light() + xlim(-2,2)

#Models per capita

models4_mayor <- (list(
  "Any Damaged Housing" = lm_robust(morena_change_mayor_15_18~damaged_houses_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(morena_change_mayor_15_18~damaged_houses_highrisk_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(morena_change_mayor_15_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('damaged_houses_1k' = 'Damage Per Capita',
               'damaged_houses_highrisk_1k' = 'Damage Per Capita',
               'damaged_multifamily_1k' = 'Damage Per Capita')

plot4_mayor <- modelplot(models4_mayor, 
                         coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                         coef_map = coefnames)+  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Challenger (MORENA)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + 
  theme_light() + xlim(-0.012, 0.012)

plot1_mayor + plot2_mayor + plot3_mayor + plot4_mayor +  plot_layout(nrow = 2)

###############################################
############ Main Text -- Figure 8 ############                     
###############################################

#Figure 8.  Relationship between damage associated with the earthquake and the electoral performance of the PRD and MORENA candidates for deputy.

#City Incumbent

#Models with distance

models1_deputies <- (list(
  "Any Damaged Housing" = lm_robust(prd_change_deputies_15_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(prd_change_deputies_15_18~distance_damaged_houses_highrisk+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(prd_change_deputies_15_18~distance_damaged_multifamily+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('distance_damaged_houses' = 'Distance to Damage',
               'distance_damaged_houses_highrisk' = 'Distance to Damage',
               'distance_damaged_multifamily' = 'Distance to Damage')

plot1_deputies <- modelplot(models1_deputies, 
                            coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                            coef_map = coefnames) +  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Incumbent (PRD)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + theme_light() + xlim(-2.5,2.5)

#Models per capita

models2_deputies <- (list(
  "Any Damaged Housing" = lm_robust(prd_change_deputies_15_18~damaged_houses_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(prd_change_deputies_15_18~damaged_houses_highrisk_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(prd_change_deputies_15_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('damaged_houses_1k' = 'Damage Per Capita',
               'damaged_houses_highrisk_1k' = 'Damage Per Capita',
               'damaged_multifamily_1k' = 'Damage Per Capita')

plot2_deputies <- modelplot(models2_deputies, 
                            coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                            coef_map = coefnames)+  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Incumbent (PRD)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + 
  theme_light() + xlim(-0.012, 0.012)

#City Challenger

#Models with distance

models3_deputies <- (list(
  "Any Damaged Housing" = lm_robust(morena_change_deputies_15_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(morena_change_deputies_15_18~distance_damaged_houses_highrisk+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(morena_change_deputies_15_18~distance_damaged_multifamily+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('distance_damaged_houses' = 'Distance to Damage',
               'distance_damaged_houses_highrisk' = 'Distance to Damage',
               'distance_damaged_multifamily' = 'Distance to Damage')

plot3_deputies <- modelplot(models3_deputies, 
                            coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                            coef_map = coefnames)+  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Challenger (MORENA)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + theme_light() + xlim(-2.5,2.5)

#Models per capita

models4_deputies <- (list(
  "Any Damaged Housing" = lm_robust(morena_change_deputies_15_18~damaged_houses_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "High Risk House" = lm_robust(morena_change_deputies_15_18~damaged_houses_highrisk_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  "Major Building" = lm_robust(morena_change_deputies_15_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('damaged_houses_1k' = 'Damage Per Capita',
               'damaged_houses_highrisk_1k' = 'Damage Per Capita',
               'damaged_multifamily_1k' = 'Damage Per Capita')

plot4_deputies <- modelplot(models4_deputies, 
                            coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|Intercept"),
                            coef_map = coefnames) +  
  geom_vline(xintercept = 0, lty="longdash") +  
  labs(title = "City Challenger (MORENA)") +
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + 
  theme_light() + xlim(-0.012, 0.012)

plot1_deputies + plot2_deputies + plot3_deputies + plot4_deputies +  plot_layout(nrow = 2)

###############################################
############ Main Text -- Figure 9 ############                     
###############################################

#Figure 9.  Relationship between damage associated with the earthquake, electoral performance of the PRD and MORENA candidates, and distribution of disaster-relief policies. 

models2_rent <- (list(
  "Mayors PRD" = lm_robust(prd_change_mayor_15_18~housingrelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Mayors MORENA" = lm_robust(morena_change_mayor_15_18~housingrelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Deputies PRD" = lm_robust(prd_change_deputies_15_18~housingrelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Deputies MORENA" = lm_robust(morena_change_deputies_15_18~housingrelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Governor PRD" = lm_robust(prd_change_gov_12_18~housingrelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('housingrelief_1k' = 'Reconstruction Credits per 1000')

plot2_rent <- modelplot(models2_rent, 
                        coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|damaged_multifamily_1k|Intercept"),
                        coef_map = coefnames) +
  labs(title = "Reconstruction Credits per 1000 Inhabitants")  +  
  geom_vline(xintercept = 0, lty="longdash") +  
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + 
  theme_light() + xlim(-0.002, 0.002)

#Procuraduria aid

models2_procu <- (list(
  "Mayors PRD" = lm_robust(prd_change_mayor_15_18~procuraduriarelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Mayors MORENA" = lm_robust(morena_change_mayor_15_18~procuraduriarelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Deputies PRD" = lm_robust(prd_change_deputies_15_18~procuraduriarelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Deputies MORENA" = lm_robust(morena_change_deputies_15_18~procuraduriarelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata"),
  
  "Governor PRD" = lm_robust(prd_change_gov_12_18~procuraduriarelief_1k+damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master, se_type = "stata")))

coefnames <- c('procuraduriarelief_1k' = 'Risk Reduction per 1000')

plot2_procu <- modelplot(models2_procu, 
                         coef_omit = c("as.factor|president2018_share_morena|seguro_popular|unemployment|literacy|urbage|damaged_multifamily_1k|Intercept"),
                         coef_map = coefnames) +
  labs(title = "Risk Reduction per 1000 Inhabitants")  +  
  geom_vline(xintercept = 0, lty="longdash") +  
  #theme(axis.text.y=element_blank(),  axis.ticks.y=element_blank()) + 
  scale_colour_grey(start = 0, end = .7) + guides(colour = guide_legend(reverse=TRUE)) + 
  theme_light() + xlim(-0.002, 0.002)

plot2_procu + plot2_rent +  plot_layout(nrow = 2)

###############################################
############ Main Text -- Table 1 #############                     
###############################################

#Table 1.  Sensitivity analysis for main results (deputies/distance-based damage).

model1 <- lm(prd_change_deputies_15_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master)

sensitivity1 <- sensemakr(model = model1, 
                          treatment = "distance_damaged_houses",
                          benchmark_covariates = "president2018_share_morena",
                          kd = 1:2)

sensemakr::ovb_minimal_reporting(sensitivity1, format = "latex")

###############################################
############ Main Text -- Table 2 #############                     
###############################################

#Table 2.  Sensitivity analysis for main results (deputies/per capita damage).

model1b <- lm(prd_change_deputies_15_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master)

sensitivity1b <- sensemakr(model = model1b, 
                           treatment = "damaged_multifamily_1k",
                           benchmark_covariates = "president2018_share_morena",
                           kd = 1:2)

sensemakr::ovb_minimal_reporting(sensitivity1b, format = "latex")

###############################################
############ Main Text -- Table 3 #############                     
###############################################

#Table 3.  Sensitivity analysis for main results (mayors/distance-based damage).

model2 <- lm(prd_change_mayor_15_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master)

sensitivity2 <- sensemakr(model = model2, 
                          treatment = "distance_damaged_houses",
                          benchmark_covariates = "president2018_share_morena",
                          kd = 1:2)

sensemakr::ovb_minimal_reporting(sensitivity2, format = "latex")

###############################################
############ Main Text -- Table 4 #############                     
###############################################

#Table 4.  Sensitivity analysis for main results (mayors/per capita damage).

model2b <- lm(prd_change_mayor_15_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master)

sensitivity2b <- sensemakr(model = model2b, 
                           treatment = "damaged_multifamily_1k",
                           benchmark_covariates = "president2018_share_morena",
                           kd = 1:2)

sensemakr::ovb_minimal_reporting(sensitivity2b, format = "latex")

###############################################
############ Main Text -- Table 5 #############                     
###############################################

#Table 5.  Sensitivity analysis for main results (governor/distance-based damage).

model3 <- lm(prd_change_gov_12_18~distance_damaged_houses+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master)

sensitivity3 <- sensemakr(model = model3, 
                          treatment = "distance_damaged_houses",
                          benchmark_covariates = "president2018_share_morena",
                          kd = 1:2)

sensemakr::ovb_minimal_reporting(sensitivity3, format = "latex")

###############################################
############ Main Text -- Table 6 #############                     
###############################################

#Table 6.  Sensitivity analysis for main results (governor/per capita damage).

model3b <- lm(prd_change_gov_12_18~damaged_multifamily_1k+as.factor(zona.geotec)+urbage+literacy_15+unemployment+seguro_popular+president2018_share_morena+as.factor(name.del), data = master)

sensitivity3b <- sensemakr(model = model3b, 
                           treatment = "damaged_multifamily_1k",
                           benchmark_covariates = "president2018_share_morena",
                           kd = 1:2)

sensemakr::ovb_minimal_reporting(sensitivity3b, format = "latex")
